# PEDL_4B.R

# Contains: Code to generate Figure 4, Table 2, and Figure A4

# Load libraries
library(pacman)
pacman::p_load(ggplot2, patchwork)

# Set working directory (YOUR FILEPATH HERE)
setwd("")

#######################################################################
#######################################################################

# Import data
primary <- read.csv(file = "primary_cleaned.csv", header = TRUE)
paircoef <- read.csv(file = "PairwiseComparisons.csv", header = TRUE)

#######################################################################
#######################################################################

# Clean arm variable
primary$arm.factor <- factor(primary$arm, levels = c("0", "125", "250", "375"))
primary$arm.dollar <- factor(ifelse(primary$arm == 0, "Control", 
                                    ifelse(primary$arm == 125, "$1.7",
                                           ifelse(primary$arm == 250, "$3.4", "$5.1"))),
                             levels = c("Control", "$1.7", "$3.4", "$5.1"))

# Create primary outcome variable
primary$intervention_weight <- primary$intervention_refills2*14.2

primary$lpg_wt_used <- primary$bline_lpg + 
  ((primary$intervention_refills2 - 1)*14.2) + (14.2 - primary$eline_lpg)

# Participants did not buy any LPG, and had a greater cylinder weight 
# at endline than baseline. These have been recoded to 0.
primary$lpg_wt_used <- ifelse(primary$lpg_wt_used < 0, 0, primary$lpg_wt_used)

# replace suspension_refills2 = 0 if suspension_refills2 == .
# gen overall_weight = 
primary$suspension_refills2 <- ifelse(is.na(primary$suspension_refills2), 0, primary$suspension_refills2)
primary$lpg_wt_used_wsus <- ((primary$intervention_refills2 + primary$suspension_refills2)*14.2) + 
  primary$bline_lpg - primary$eline_lpg



# Figure 4 (Note: these values come from the Stata code)
mod_coef <- data.frame(
  Comp = factor(c("$1.7 - Control", "$3.4 - Control", 
                  "$5.1 - Control", "$3.4 - $1.7", 
                  "$5.1 - $1.7", "$5.1 - $3.4"), 
                levels = c("$5.1 - $3.4", "$5.1 - $1.7", 
                           "$3.4 - $1.7", "$5.1 - Control", 
                           "$3.4 - Control", "$1.7 - Control")),
  coef = paircoef$b,
  lci = paircoef$lci,
  uci = paircoef$uci
)

f4 <- ggplot(data = mod_coef, 
             aes(x=Comp, y=coef)) + 
  geom_point(position=position_dodge(0.2)) + 
  geom_errorbar(aes(ymin=lci, ymax=uci), width=.1, position=position_dodge(0.2)) + 
  scale_color_manual(values = "black") +
  geom_hline(yintercept=0, linetype=2) +
  xlab('Pairwise comparison') + 
  ylab("Coefficient (95% CI)") +
  ylim(c(-3, 22)) +
  theme_bw() +
  coord_flip() 

f4





# Table 2
aggregate( lpg_wt_used ~ arm.factor, primary, mean )
aggregate( lpg_wt_used ~ arm.factor, primary, sd )

aggregate( preEstlpgPer30_1 ~ arm.factor, primary, mean )
aggregate( preEstlpgPer30_1 ~ arm.factor, primary, sd )

aggregate( postEstlpgPer30_1 ~ arm.factor, primary, mean )
aggregate( postEstlpgPer30_1 ~ arm.factor, primary, sd )

aggregate( lpg_wt_used_wsus ~ arm.factor, primary, mean )
aggregate( lpg_wt_used_wsus ~ arm.factor, primary, sd )






# Models for pre-suspension, post-suspension, wt used w/ suspension refills
mod1 <- aov(preEstlpgPer30_1 ~ arm.factor, data = primary)
summary(mod1)
summary.lm(mod1)
mod1_Tuk <- TukeyHSD(mod1)

mod1_coef <- data.frame(
  Comp = factor(c("$1.7 - Control", "$3.4 - Control", 
                  "$5.1 - Control", "$3.4 - $1.7", 
                  "$5.1 - $1.7", "$5.1 - $3.4"), 
                levels = c("$5.1 - $3.4", "$5.1 - $1.7", 
                           "$3.4 - $1.7", "$5.1 - Control", 
                           "$3.4 - Control", "$1.7 - Control")),
  coef = mod1_Tuk$arm.factor[1:6],
  lci = mod1_Tuk$arm.factor[7:12],
  uci = mod1_Tuk$arm.factor[13:18]
)



mod2 <- aov(postEstlpgPer30_1 ~ arm.factor, data = primary)
summary(mod2)
summary.lm(mod2)
mod2_Tuk <- TukeyHSD(mod2)

mod2_coef <- data.frame(
  Comp = factor(c("$1.7 - Control", "$3.4 - Control", 
                  "$5.1 - Control", "$3.4 - $1.7", 
                  "$5.1 - $1.7", "$5.1 - $3.4"), 
                levels = c("$5.1 - $3.4", "$5.1 - $1.7", 
                           "$3.4 - $1.7", "$5.1 - Control", 
                           "$3.4 - Control", "$1.7 - Control")),
  coef = mod2_Tuk$arm.factor[1:6],
  lci = mod2_Tuk$arm.factor[7:12],
  uci = mod2_Tuk$arm.factor[13:18]
)



mod3 <- aov(lpg_wt_used_wsus ~ arm.factor, data = primary)
summary(mod3)
summary.lm(mod3)
mod3_Tuk <- TukeyHSD(mod3)



# Figure A4
fa4_1 <- ggplot(data = mod1_coef, 
              aes(x=Comp, y=coef)) + 
  geom_point(position=position_dodge(0.2)) + 
  geom_errorbar(aes(ymin=lci, ymax=uci), width=.1, position=position_dodge(0.2)) + 
  scale_color_manual(values = "black") +
  geom_hline(yintercept=0, linetype=2) +
  xlab('Pairwise comparison') + 
  ylab("Coefficient (95% CI)") +
  ylim(c(-1, 3)) +
  theme_bw() +
  coord_flip() +
  ggtitle("A. Pre-suspension weight of LPG consumed \n per 30 days")

fa4_1


fa4_2 <- ggplot(data = mod2_coef, 
              aes(x=Comp, y=coef)) + 
  geom_point(position=position_dodge(0.2)) + 
  geom_errorbar(aes(ymin=lci, ymax=uci), width=.1, position=position_dodge(0.2)) + 
  scale_color_manual(values = "black") +
  geom_hline(yintercept=0, linetype=2) +
  xlab('Pairwise comparison') + 
  ylab("Coefficient (95% CI)") +
  ylim(c(-1, 3)) +
  theme_bw() +
  coord_flip() +
  ggtitle("B. Post-suspension weight of LPG consumed \n per 30 days")

fa4_2

fa4 <- fa4_1 + fa4_2
fa4


